﻿/*
 Copyright (c) 2011 Seth Juarez

 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights
 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 copies of the Software, and to permit persons to whom the Software is
 furnished to do so, subject to the following conditions:

 The above copyright notice and this permission notice shall be included in
 all copies or substantial portions of the Software.

 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
*/

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ml.Math
{
    /// <summary>
    /// Eigenvalues and eigenvectors of a real matrix. 
    ///
    /// If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
    /// diagonal and the eigenvector matrix V is orthogonal.
    /// I.e. A = V.times(D.times(V.transpose())) and 
    /// V.times(V.transpose()) equals the identity matrix.
    ///
    /// If A is not symmetric, then the eigenvalue matrix D is block diagonal
    /// with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
    /// lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda].  The
    /// columns of V represent the eigenvectors in the sense that A*V = V*D,
    /// i.e. A.times(V) equals V.times(D).  The matrix V may be badly
    /// conditioned, or even singular, so the validity of the equation
    /// A = V*D*inverse(V) depends upon V.cond().
    /// 
    /// This class was ported from Java to C# from JAMA: http://math.nist.gov/javanumerics/jama/
    /// Original Source File: EigenvalueDecomposition.java
    /// </summary>
    internal class EigenDecomposition
    {

        #region Class variables

        /// <summary>
        /// Row and column dimension (square matrix).
        /// </summary>
        private int n { get; set; }

        /// <summary>
        /// internal symmetry flag.
        /// </summary>
        private bool sym { get; set; }

        /// <summary>
        /// Imaginary eigenvalues
        /// </summary>
        public Vector e { get; private set; }

        /// <summary>
        /// Real eigenvalues
        /// </summary>
        public Vector d { get; private set; }

        /// <summary>
        /// Block diagonal eigenvalue matrix
        /// (Computed on Demand)
        /// </summary>
        private Matrix _D;

        /// <summary>
        /// Block diagonal eigenvalue matrix
        /// </summary>
        public Matrix D
        {
            get
            {
                if (_D == null)
                {
                    _D = Matrix.Zeros(n);
                    for (int i = 0; i < n; i++)
                    {
                        _D[i, i] = d[i];
                        if (e[i] > 0)
                            _D[i, i + 1] = e[i];
                        else if (e[i] < 0)
                            _D[i, i - 1] = e[i];
                    }
                }

                return _D;
            }
        }

        /// <summary>
        /// Eigenvector Matrix
        /// </summary>
        public Matrix V { get; private set; }

        /// <summary>
        /// Array for internal storage of nonsymmetric Hessenberg form.
        /// </summary>
        private Matrix H { get; set; }

        /// <summary>
        /// Working storage for nonsymmetric algorithm.
        /// </summary>
        private Vector ort { get; set; }

        #endregion

        #region Private Methods
        /// <summary>
        /// Symmetric Householder reduction to tridiagonal form.
        /// </summary>
        private void tred2()
        {

            //  This is derived from the Algol procedures tred2 by
            //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
            //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutine in EISPACK.

            for (int j = 0; j < n; j++)
                d[j] = V[n - 1, j];

            // Householder reduction to tridiagonal form.

            for (int i = n - 1; i > 0; i--)
            {

                // Scale to avoid under/overflow.

                double scale = 0.0;
                double h = 0.0;
                for (int k = 0; k < i; k++)
                    scale = scale + System.Math.Abs(d[k]);

                if (scale == 0.0)
                {
                    e[i] = d[i - 1];
                    for (int j = 0; j < i; j++)
                    {
                        d[j] = V[i - 1, j];
                        V[i, j] = 0.0;
                        V[j, i] = 0.0;
                    }
                }
                else
                {

                    // Generate Householder vector.
                    for (int k = 0; k < i; k++)
                    {
                        d[k] /= scale;
                        h += d[k] * d[k];
                    }

                    double f = d[i - 1];
                    double g = System.Math.Sqrt(h);

                    if (f > 0)
                        g = -g;

                    e[i] = scale * g;
                    h = h - f * g;
                    d[i - 1] = f - g;

                    for (int j = 0; j < i; j++)
                        e[j] = 0.0;

                    // Apply similarity transformation to remaining columns.

                    for (int j = 0; j < i; j++)
                    {
                        f = d[j];
                        V[j, i] = f;
                        g = e[j] + V[j, j] * f;
                        for (int k = j + 1; k <= i - 1; k++)
                        {
                            g += V[k, j] * d[k];
                            e[k] += V[k, j] * f;
                        }
                        e[j] = g;
                    }

                    f = 0.0;
                    for (int j = 0; j < i; j++)
                    {
                        e[j] /= h;
                        f += e[j] * d[j];
                    }
                    double hh = f / (h + h);
                    for (int j = 0; j < i; j++)
                        e[j] -= hh * d[j];

                    for (int j = 0; j < i; j++)
                    {
                        f = d[j];
                        g = e[j];

                        for (int k = j; k <= i - 1; k++)
                            V[k, j] -= (f * e[k] + g * d[k]);

                        d[j] = V[i - 1, j];
                        V[i, j] = 0.0;
                    }
                }
                d[i] = h;
            }

            // Accumulate transformations.
            for (int i = 0; i < n - 1; i++)
            {
                V[n - 1, i] = V[i, i];
                V[i, i] = 1.0;
                double h = d[i + 1];
                if (h != 0.0)
                {
                    for (int k = 0; k <= i; k++)
                        d[k] = V[k, i + 1] / h;

                    for (int j = 0; j <= i; j++)
                    {
                        double g = 0.0;
                        for (int k = 0; k <= i; k++)
                            g += V[k, i + 1] * V[k, j];

                        for (int k = 0; k <= i; k++)
                            V[k, j] -= g * d[k];

                    }
                }

                for (int k = 0; k <= i; k++)
                    V[k, i + 1] = 0.0;

            }
            for (int j = 0; j < n; j++)
            {
                d[j] = V[n - 1, j];
                V[n - 1, j] = 0.0;
            }
            V[n - 1, n - 1] = 1.0;
            e[0] = 0.0;
        }

        /// <summary>
        /// Symmetric tridiagonal QL algorithm.
        /// </summary>
        private void tql2()
        {

            //  This is derived from the Algol procedures tql2, by
            //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
            //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutine in EISPACK.

            for (int i = 1; i < n; i++)
                e[i - 1] = e[i];

            e[n - 1] = 0.0;

            double f = 0.0;
            double tst1 = 0.0;
            double eps = System.Math.Pow(2.0, -52.0);
            for (int l = 0; l < n; l++)
            {

                // Find small subdiagonal element

                tst1 = System.Math.Max(tst1, System.Math.Abs(d[l]) + System.Math.Abs(e[l]));
                int m = l;
                while (m < n)
                {
                    if (System.Math.Abs(e[m]) <= eps * tst1)
                        break;

                    m++;
                }

                // If m == l, d[l] is an eigenvalue,
                // otherwise, iterate.

                if (m > l)
                {
                    int iter = 0;
                    do
                    {
                        iter = iter + 1;  // (Could check iteration count here.)

                        // Compute implicit shift

                        double g = d[l];
                        double p = (d[l + 1] - g) / (2.0 * e[l]);
                        double r = new Vector(new double[] { p, 1.0 }).Norm(); //Helper.Hypot(p, 1.0);
                        if (p < 0)
                            r = -r;

                        d[l] = e[l] / (p + r);
                        d[l + 1] = e[l] * (p + r);
                        double dl1 = d[l + 1];
                        double h = g - d[l];
                        for (int i = l + 2; i < n; i++)
                            d[i] -= h;

                        f = f + h;

                        // Implicit QL transformation.
                        p = d[m];
                        double c = 1.0;
                        double c2 = c;
                        double c3 = c;
                        double el1 = e[l + 1];
                        double s = 0.0;
                        double s2 = 0.0;
                        for (int i = m - 1; i >= l; i--)
                        {
                            c3 = c2;
                            c2 = c;
                            s2 = s;
                            g = c * e[i];
                            h = c * p;
                            r = new Vector(new double[] { p, e[i] }).Norm(); //Helper.Hypot(p, e[i]);
                            e[i + 1] = s * r;
                            s = e[i] / r;
                            c = p / r;
                            p = c * d[i] - s * g;
                            d[i + 1] = h + s * (c * g + s * d[i]);

                            // Accumulate transformation.
                            for (int k = 0; k < n; k++)
                            {
                                h = V[k, i + 1];
                                V[k, i + 1] = s * V[k, i] + c * h;
                                V[k, i] = c * V[k, i] - s * h;
                            }
                        }

                        p = -s * s2 * c3 * el1 * e[l] / dl1;
                        e[l] = s * p;
                        d[l] = c * p;

                        // Check for convergence.

                    } while (System.Math.Abs(e[l]) > eps * tst1);
                }

                d[l] = d[l] + f;
                e[l] = 0.0;
            }

            // Sort eigenvalues and corresponding vectors.

            for (int i = 0; i < n - 1; i++)
            {
                int k = i;
                double p = d[i];
                for (int j = i + 1; j < n; j++)
                {
                    if (d[j] < p)
                    {
                        k = j;
                        p = d[j];
                    }
                }
                if (k != i)
                {
                    d[k] = d[i];
                    d[i] = p;
                    for (int j = 0; j < n; j++)
                    {
                        p = V[j, i];
                        V[j, i] = V[j, k];
                        V[j, k] = p;
                    }
                }
            }
        }

        /// <summary>
        /// Nonsymmetric reduction to Hessenberg form.
        /// </summary>
        private void orthes()
        {

            //  This is derived from the Algol procedures orthes and ortran,
            //  by Martin and Wilkinson, Handbook for Auto. Comp.,
            //  Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutines in EISPACK.

            int low = 0;
            int high = n - 1;

            for (int m = low + 1; m <= high - 1; m++)
            {

                // Scale column.

                double scale = 0.0;
                for (int i = m; i <= high; i++)
                {
                    scale = scale + System.Math.Abs(H[i, m - 1]);
                }
                if (scale != 0.0)
                {

                    // Compute Householder transformation.

                    double h = 0.0;
                    for (int i = high; i >= m; i--)
                    {
                        ort[i] = H[i, m - 1] / scale;
                        h += ort[i] * ort[i];
                    }
                    double g = System.Math.Sqrt(h);
                    if (ort[m] > 0)
                    {
                        g = -g;
                    }
                    h = h - ort[m] * g;
                    ort[m] = ort[m] - g;

                    // Apply Householder similarity transformation
                    // H = (I-u*u'/h)*H*(I-u*u')/h)

                    for (int j = m; j < n; j++)
                    {
                        double f = 0.0;
                        for (int i = high; i >= m; i--)
                        {
                            f += ort[i] * H[i, j];
                        }
                        f = f / h;
                        for (int i = m; i <= high; i++)
                        {
                            H[i, j] -= f * ort[i];
                        }
                    }

                    for (int i = 0; i <= high; i++)
                    {
                        double f = 0.0;
                        for (int j = high; j >= m; j--)
                        {
                            f += ort[j] * H[i, j];
                        }
                        f = f / h;
                        for (int j = m; j <= high; j++)
                        {
                            H[i, j] -= f * ort[j];
                        }
                    }
                    ort[m] = scale * ort[m];
                    H[m, m - 1] = scale * g;
                }
            }

            // Accumulate transformations (Algol's ortran).

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    V[i, j] = (i == j ? 1.0 : 0.0);
                }
            }

            for (int m = high - 1; m >= low + 1; m--)
            {
                if (H[m, m - 1] != 0.0)
                {
                    for (int i = m + 1; i <= high; i++)
                    {
                        ort[i] = H[i, m - 1];
                    }
                    for (int j = m; j <= high; j++)
                    {
                        double g = 0.0;
                        for (int i = m; i <= high; i++)
                        {
                            g += ort[i] * V[i, j];
                        }
                        // Double division avoids possible underflow
                        g = (g / ort[m]) / H[m, m - 1];
                        for (int i = m; i <= high; i++)
                        {
                            V[i, j] += g * ort[i];
                        }
                    }
                }
            }
        }

        // Complex scalar division.
        private double cdivr, cdivi;
        private void cdiv(double xr, double xi, double yr, double yi)
        {
            double r, d;
            if (System.Math.Abs(yr) > System.Math.Abs(yi))
            {
                r = yi / yr;
                d = yr + r * yi;
                cdivr = (xr + r * xi) / d;
                cdivi = (xi - r * xr) / d;
            }
            else
            {
                r = yr / yi;
                d = yi + r * yr;
                cdivr = (r * xr + xi) / d;
                cdivi = (r * xi - xr) / d;
            }
        }

        /// <summary>
        /// Nonsymmetric reduction from Hessenberg to real Schur form.
        /// </summary>
        private void hqr2()
        {

            //  This is derived from the Algol procedure hqr2,
            //  by Martin and Wilkinson, Handbook for Auto. Comp.,
            //  Vol.ii-Linear Algebra, and the corresponding
            //  Fortran subroutine in EISPACK.

            // Initialize

            int nn = this.n;
            int n = nn - 1;
            int low = 0;
            int high = nn - 1;
            double eps = System.Math.Pow(2.0, -52.0);
            double exshift = 0.0;
            double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;

            // Store roots isolated by balanc and compute matrix norm

            double norm = 0.0;
            for (int i = 0; i < nn; i++)
            {
                if (i < low | i > high)
                {
                    d[i] = H[i, i];
                    e[i] = 0.0;
                }
                for (int j = System.Math.Max(i - 1, 0); j < nn; j++)
                {
                    norm = norm + System.Math.Abs(H[i, j]);
                }
            }

            // Outer loop over eigenvalue index

            int iter = 0;
            while (n >= low)
            {

                // Look for single small sub-diagonal element

                int l = n;
                while (l > low)
                {
                    s = System.Math.Abs(H[l - 1, l - 1]) + System.Math.Abs(H[l, l]);
                    if (s == 0.0)
                        s = norm;
                    if (System.Math.Abs(H[l, l - 1]) < eps * s)
                        break;
                    l--;
                }

                // Check for convergence
                // One root found

                if (l == n)
                {
                    H[n, n] = H[n, n] + exshift;
                    d[n] = H[n, n];
                    e[n] = 0.0;
                    n--;
                    iter = 0;

                    // Two roots found

                }
                else if (l == n - 1)
                {
                    w = H[n, n - 1] * H[n - 1, n];
                    p = (H[n - 1, n - 1] - H[n, n]) / 2.0;
                    q = p * p + w;
                    z = System.Math.Sqrt(System.Math.Abs(q));
                    H[n, n] = H[n, n] + exshift;
                    H[n - 1, n - 1] = H[n - 1, n - 1] + exshift;
                    x = H[n, n];

                    // Real pair

                    if (q >= 0)
                    {
                        if (p >= 0)
                            z = p + z;
                        else
                            z = p - z;
                        d[n - 1] = x + z;
                        d[n] = d[n - 1];
                        if (z != 0.0)
                            d[n] = x - w / z;
                        e[n - 1] = 0.0;
                        e[n] = 0.0;
                        x = H[n, n - 1];
                        s = System.Math.Abs(x) + System.Math.Abs(z);
                        p = x / s;
                        q = z / s;
                        r = System.Math.Sqrt(p * p + q * q);
                        p = p / r;
                        q = q / r;

                        // Row modification
                        for (int j = n - 1; j < nn; j++)
                        {
                            z = H[n - 1, j];
                            H[n - 1, j] = q * z + p * H[n, j];
                            H[n, j] = q * H[n, j] - p * z;
                        }

                        // Column modification
                        for (int i = 0; i <= n; i++)
                        {
                            z = H[i, n - 1];
                            H[i, n - 1] = q * z + p * H[i, n];
                            H[i, n] = q * H[i, n] - p * z;
                        }

                        // Accumulate transformations
                        for (int i = low; i <= high; i++)
                        {
                            z = V[i, n - 1];
                            V[i, n - 1] = q * z + p * V[i, n];
                            V[i, n] = q * V[i, n] - p * z;
                        }

                        // Complex pair

                    }
                    else
                    {
                        d[n - 1] = x + p;
                        d[n] = x + p;
                        e[n - 1] = z;
                        e[n] = -z;
                    }
                    n = n - 2;
                    iter = 0;

                    // No convergence yet

                }
                else
                {

                    // Form shift
                    x = H[n, n];
                    y = 0.0;
                    w = 0.0;
                    if (l < n)
                    {
                        y = H[n - 1, n - 1];
                        w = H[n, n - 1] * H[n - 1, n];
                    }

                    // Wilkinson's original ad hoc shift
                    if (iter == 10)
                    {
                        exshift += x;
                        for (int i = low; i <= n; i++)
                            H[i, i] -= x;
                        s = System.Math.Abs(H[n, n - 1]) + System.Math.Abs(H[n - 1, n - 2]);
                        x = y = 0.75 * s;
                        w = -0.4375 * s * s;
                    }

                    // MATLAB's new ad hoc shift
                    if (iter == 30)
                    {
                        s = (y - x) / 2.0;
                        s = s * s + w;
                        if (s > 0)
                        {
                            s = System.Math.Sqrt(s);

                            if (y < x)
                                s = -s;

                            s = x - w / ((y - x) / 2.0 + s);

                            for (int i = low; i <= n; i++)
                                H[i, i] -= s;

                            exshift += s;
                            x = y = w = 0.964;
                        }
                    }

                    iter = iter + 1;   // (Could check iteration count here.)

                    // Look for two consecutive small sub-diagonal elements
                    int m = n - 2;
                    while (m >= l)
                    {
                        z = H[m, m];
                        r = x - z;
                        s = y - z;
                        p = (r * s - w) / H[m + 1, m] + H[m, m + 1];
                        q = H[m + 1, m + 1] - z - r - s;
                        r = H[m + 2, m + 1];
                        s = System.Math.Abs(p) + System.Math.Abs(q) + System.Math.Abs(r);
                        p = p / s;
                        q = q / s;
                        r = r / s;
                        if (m == l)
                            break;

                        if (System.Math.Abs(H[m, m - 1]) * (System.Math.Abs(q) + System.Math.Abs(r)) <
                           eps * (System.Math.Abs(p) * (System.Math.Abs(H[m - 1, m - 1]) + System.Math.Abs(z) +
                           System.Math.Abs(H[m + 1, m + 1])))) break;

                        m--;
                    }

                    for (int i = m + 2; i <= n; i++)
                    {
                        H[i, i - 2] = 0.0;
                        if (i > m + 2)
                            H[i, i - 3] = 0.0;
                    }

                    // Double QR step involving rows l:n and columns m:n
                    for (int k = m; k <= n - 1; k++)
                    {
                        bool notlast = (k != n - 1);
                        if (k != m)
                        {
                            p = H[k, k - 1];
                            q = H[k + 1, k - 1];
                            r = (notlast ? H[k + 2, k - 1] : 0.0);
                            x = System.Math.Abs(p) + System.Math.Abs(q) + System.Math.Abs(r);
                            if (x != 0.0)
                            {
                                p = p / x;
                                q = q / x;
                                r = r / x;
                            }
                        }
                        if (x == 0.0)
                            break;

                        s = System.Math.Sqrt(p * p + q * q + r * r);
                        if (p < 0)
                            s = -s;

                        if (s != 0)
                        {
                            if (k != m)
                                H[k, k - 1] = -s * x;
                            else if (l != m)
                                H[k, k - 1] = -H[k, k - 1];

                            p = p + s;
                            x = p / s;
                            y = q / s;
                            z = r / s;
                            q = q / p;
                            r = r / p;

                            // Row modification
                            for (int j = k; j < nn; j++)
                            {
                                p = H[k, j] + q * H[k + 1, j];
                                if (notlast)
                                {
                                    p = p + r * H[k + 2, j];
                                    H[k + 2, j] = H[k + 2, j] - p * z;
                                }
                                H[k, j] = H[k, j] - p * x;
                                H[k + 1, j] = H[k + 1, j] - p * y;
                            }

                            // Column modification
                            for (int i = 0; i <= System.Math.Min(n, k + 3); i++)
                            {
                                p = x * H[i, k] + y * H[i, k + 1];
                                if (notlast)
                                {
                                    p = p + z * H[i, k + 2];
                                    H[i, k + 2] = H[i, k + 2] - p * r;
                                }
                                H[i, k] = H[i, k] - p;
                                H[i, k + 1] = H[i, k + 1] - p * q;
                            }

                            // Accumulate transformations
                            for (int i = low; i <= high; i++)
                            {
                                p = x * V[i, k] + y * V[i, k + 1];
                                if (notlast)
                                {
                                    p = p + z * V[i, k + 2];
                                    V[i, k + 2] = V[i, k + 2] - p * r;
                                }
                                V[i, k] = V[i, k] - p;
                                V[i, k + 1] = V[i, k + 1] - p * q;
                            }
                        }  // (s != 0)
                    }  // k loop
                }  // check convergence
            }  // while (n >= low)

            // Backsubstitute to find vectors of upper triangular form
            if (norm == 0.0)
                return;

            for (n = nn - 1; n >= 0; n--)
            {
                p = d[n];
                q = e[n];

                // Real vector
                if (q == 0)
                {
                    int l = n;
                    H[n, n] = 1.0;
                    for (int i = n - 1; i >= 0; i--)
                    {
                        w = H[i, i] - p;
                        r = 0.0;
                        for (int j = l; j <= n; j++)
                            r = r + H[i, j] * H[j, n];

                        if (e[i] < 0.0)
                        {
                            z = w;
                            s = r;
                        }
                        else
                        {
                            l = i;
                            if (e[i] == 0.0)
                            {
                                if (w != 0.0)
                                    H[i, n] = -r / w;
                                else
                                    H[i, n] = -r / (eps * norm);

                                // Solve real equations
                            }
                            else
                            {
                                x = H[i, i + 1];
                                y = H[i + 1, i];
                                q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
                                t = (x * s - z * r) / q;
                                H[i, n] = t;
                                if (System.Math.Abs(x) > System.Math.Abs(z))
                                    H[i + 1, n] = (-r - w * t) / x;
                                else
                                    H[i + 1, n] = (-s - y * t) / z;
                            }

                            // Overflow control

                            t = System.Math.Abs(H[i, n]);
                            if ((eps * t) * t > 1)
                                for (int j = i; j <= n; j++)
                                    H[j, n] = H[j, n] / t;
                        }
                    }

                    // Complex vector

                }
                else if (q < 0)
                {
                    int l = n - 1;

                    // Last vector component imaginary so matrix is triangular

                    if (System.Math.Abs(H[n, n - 1]) > System.Math.Abs(H[n - 1, n]))
                    {
                        H[n - 1, n - 1] = q / H[n, n - 1];
                        H[n - 1, n] = -(H[n, n] - p) / H[n, n - 1];
                    }
                    else
                    {
                        cdiv(0.0, -H[n - 1, n], H[n - 1, n - 1] - p, q);
                        H[n - 1, n - 1] = cdivr;
                        H[n - 1, n] = cdivi;
                    }
                    H[n, n - 1] = 0.0;
                    H[n, n] = 1.0;
                    for (int i = n - 2; i >= 0; i--)
                    {
                        double ra, sa, vr, vi;
                        ra = 0.0;
                        sa = 0.0;
                        for (int j = l; j <= n; j++)
                        {
                            ra = ra + H[i, j] * H[j, n - 1];
                            sa = sa + H[i, j] * H[j, n];
                        }
                        w = H[i, i] - p;

                        if (e[i] < 0.0)
                        {
                            z = w;
                            r = ra;
                            s = sa;
                        }
                        else
                        {
                            l = i;
                            if (e[i] == 0)
                            {
                                cdiv(-ra, -sa, w, q);
                                H[i, n - 1] = cdivr;
                                H[i, n] = cdivi;
                            }
                            else
                            {
                                // Solve complex equations
                                x = H[i, i + 1];
                                y = H[i + 1, i];
                                vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
                                vi = (d[i] - p) * 2.0 * q;
                                if (vr == 0.0 & vi == 0.0)
                                {
                                    vr = eps * norm * (System.Math.Abs(w) + System.Math.Abs(q) +
                                    System.Math.Abs(x) + System.Math.Abs(y) + System.Math.Abs(z));
                                }
                                cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
                                H[i, n - 1] = cdivr;
                                H[i, n] = cdivi;
                                if (System.Math.Abs(x) > (System.Math.Abs(z) + System.Math.Abs(q)))
                                {
                                    H[i + 1, n - 1] = (-ra - w * H[i, n - 1] + q * H[i, n]) / x;
                                    H[i + 1, n] = (-sa - w * H[i, n] - q * H[i, n - 1]) / x;
                                }
                                else
                                {
                                    cdiv(-r - y * H[i, n - 1], -s - y * H[i, n], z, q);
                                    H[i + 1, n - 1] = cdivr;
                                    H[i + 1, n] = cdivi;
                                }
                            }

                            // Overflow control
                            t = System.Math.Max(System.Math.Abs(H[i, n - 1]), System.Math.Abs(H[i, n]));
                            if ((eps * t) * t > 1)
                            {
                                for (int j = i; j <= n; j++)
                                {
                                    H[j, n - 1] = H[j, n - 1] / t;
                                    H[j, n] = H[j, n] / t;
                                }
                            }
                        }
                    }
                }
            }

            // Vectors of isolated roots
            for (int i = 0; i < nn; i++)
                if (i < low | i > high)
                    for (int j = i; j < nn; j++)
                        V[i, j] = H[i, j];

            // Back transformation to get eigenvectors of original matrix
            for (int j = nn - 1; j >= low; j--)
            {
                for (int i = low; i <= high; i++)
                {
                    z = 0.0;
                    for (int k = low; k <= System.Math.Min(j, high); k++)
                        z = z + V[i, k] * H[k, j];
                    V[i, j] = z;
                }
            }
        }
        #endregion

        /// <summary>
        /// .ctor
        /// Check for symmetry, then construct the eigenvalue decomposition.
        /// Solutions are stored in D, V
        /// </summary>
        public EigenDecomposition(Matrix A)
        {
            n = A.Cols;
            V = Matrix.Zeros(n);
            d = Vector.Zeros(n);
            e = Vector.Zeros(n);

            sym = true;
            for (int j = 0; (j < n) & sym; j++)
                for (int i = 0; (i < n) & sym; i++)
                    sym = (A[i, j] == A[j, i]);

            if (sym)
            {
                for (int i = 0; i < n; i++)
                    for (int j = 0; j < n; j++)
                        V[i, j] = A[i, j];

                // Tridiagonalize.
                tred2();

                // Diagonalize.
                tql2();

            }
            else
            {
                H = Matrix.Zeros(n);
                ort = Vector.Zeros(n);

                for (int j = 0; j < n; j++)
                    for (int i = 0; i < n; i++)
                        H[i, j] = A[i, j];

                // Reduce to Hessenberg form.
                orthes();

                // Reduce Hessenberg to real Schur form.
                hqr2();
            }
        }
        
    }
}
